﻿#pragma once

#include <xy/geom/bezier.h>

namespace xy
{

namespace geom
{

template <std::size_t N, bool Rat = false> class bezier_surface : public bezier<2, N, Rat>
{
    using vecd_t = vec<2, double>;       // 参数域向量
    using vecr_t = vec<N + Rat, double>; // (wP,w)
    using vec_t = vec<N, double>;        // P
    using bezier<2, N, Rat>::m_cpts;

  public:
    using bezier<2, N, Rat>::degree;
    using bezier<2, N, Rat>::set_poles;
    using bezier<2, N, Rat>::pole_num;
    using bezier<2, N, Rat>::weight;
    using bezier<2, N, Rat>::set_weight;
    using bezier<2, N, Rat>::bound;

    /**
     * @brief BezierSurface 构造函数
     *
     * @param[in] cpts
     * @param[in] uDeg
     * @param[in] vDeg
     * @param[in] normalize
     */
    bezier_surface(std::vector<vecr_t> &&cpts, std::size_t uDeg, std::size_t vDeg, bool normalize = true)
        : bezier<2, N, Rat>(math::ndarray<2, vecr_t>{{std::move(cpts)}}, {{uDeg, vDeg}}, normalize)
    {
    }

    /**
     * @brief 计算曲线上点的位置
     *
     * @param[in] param
     * @return vec_t
     */
    vec_t d0(const vecd_t &param) const override
    {
        auto u = param[0], v = param[1];
        auto p = de_casteljau2(m_cpts.data(), pole_num(0), pole_num(1), u, v);
        if constexpr (Rat)
            return homogeneous_to_cartesian(p);
        else
            return p;
    }

    /**
     * @brief 计算曲线上点的位置和导数
     *
     * @param[in] param
     * @return std::array<vec_t, 3>
     */
    std::array<vec_t, 3> d1(const vecd_t &param) const override
    {
        std::size_t ud = degree(0);
        std::size_t vd = degree(1);
        std::vector<double> uKnots((ud + 1) * 2);
        std::vector<double> vKnots((vd + 1) * 2);
        std::fill(uKnots.begin(), uKnots.begin() + ud + 1, 0.0);
        std::fill(uKnots.begin() + ud + 1, uKnots.end(), 1.0);
        std::fill(vKnots.begin(), vKnots.begin() + vd + 1, 0.0);
        std::fill(vKnots.begin() + vd + 1, vKnots.end(), 1.0);

        auto u = param[0], v = param[1];
        if constexpr (Rat)
        {
            auto ders = rat_surface_deriv(m_cpts.data(), uKnots.data(), vKnots.data(), uKnots.size(), vKnots.size(),
                                          degree(0), degree(1), 1, u, v);
            return {ders[0][0], ders[1][0], ders[0][1]};
        }
        else
        {
            auto ders = surface_deriv_alg1(m_cpts.data(), uKnots.data(), vKnots.data(), uKnots.size(), vKnots.size(),
                                           degree(0), degree(1), 1, u, v);
            return {ders[0][0], ders[1][0], ders[0][1]};
        }
    }

    /**
     * @brief 计算曲线上点的位置、导数和二阶导数
     *
     * @param[in] param
     * @return std::array<vec_t, 6>
     */
    std::array<vec_t, 6> d2(const vecd_t &param) const override
    {
        std::size_t ud = degree(0);
        std::size_t vd = degree(1);
        std::vector<double> uKnots((ud + 1) * 2);
        std::vector<double> vKnots((vd + 1) * 2);
        std::fill(uKnots.begin(), uKnots.begin() + ud + 1, 0.0);
        std::fill(uKnots.begin() + ud + 1, uKnots.end(), 1.0);
        std::fill(vKnots.begin(), vKnots.begin() + vd + 1, 0.0);
        std::fill(vKnots.begin() + vd + 1, vKnots.end(), 1.0);

        auto u = param[0], v = param[1];
        if constexpr (Rat)
        {
            auto ders = rat_surface_deriv(m_cpts.data(), uKnots.data(), vKnots.data(), uKnots.size(), vKnots.size(),
                                          degree(0), degree(1), 2, u, v);
            return {ders[0][0], ders[1][0], ders[0][1], ders[2][0], ders[0][2], ders[1][1]};
        }
        else
        {
            auto ders = surface_deriv_alg1(m_cpts.data(), uKnots.data(), vKnots.data(), uKnots.size(), vKnots.size(),
                                           degree(0), degree(1), 2, u, v);
            return {ders[0][0], ders[1][0], ders[0][1], ders[2][0], ders[0][2], ders[1][1]};
        }
    }
};

} // namespace geom

} // namespace xy